function plt (L)

   ## usage: plt (L)
   ## 
   ## plot results

   if (!L), return ; endif

   f = dbstack ; printf("\n\t###   %s   ###\n\n", f(1).name) ;

   global PAR PRJ GVAR ADIR SDIR CDIR BATCH Judf

   Lptr = !true ; wLptr = false ;

   mth = {"J" "F" "M" "A" "M" "J" "J" "A" "S" "O" "N" "D"} ;
   P0 = 86400 ; T0 = 273.16 ; per = [2002 1 1 ; 2002 12 31] ;

   ldata = fullfile("data", PRJ, PAR.pdd) ;
   gdata = fullfile("data", PRJ, "ptr") ;
   pmkdir(pdir = fullfile(SDIR, PAR.pdd, "plt")) ;
   METAFILE = fullfile(ldata, "meta.txt") ;

   if strcmp(PAR.sim{:}), exit ; endif

   annarg = {} ; Sfun = {"mean" "max"} ;
   # annarg = {"f"} ; sfun = @(x,p) P_stat(x,p) ;

   ## plot parameters
   if BATCH
      graphics_toolkit gnuplot ; vis = "off" ;
      lw = 6*[2 1 1 1] ;
   else
      graphics_toolkit fltk ; vis = "on" ;
      lw = 1*[2 1 1 1] ;
      set(0, "defaultfigureposition", [0 1 0.01 0.01])
   endif
   pg = [2 2] ; q = 0.1 ; legpos = "northwest" ;
   cmap = [0 0 0; 0 0 1; rainbow(16)([10 2],:)] ;
   set(0, "defaultaxesfontname", "helvetica", "defaultaxesfontsize", 8, "defaultaxesfontweight", "bold") ;
   set(0, "defaulttextfontname", "helvetica", "defaulttextfontsize", 9, "defaulttextfontweight", "bold") ;

   pkg load econometrics general miscellaneous nan statistics io

   [ids, vars, names, lons, lats, alts, undefs, J] = meta(METAFILE) ;
   VAR = sunique(vars) ; VARNAME = sunique(names) ;

   if Lptr
      printf("plt: <-- %s\n", file = fullfile(ADIR, "ptrz.ob")) ;
      load(file) ;
      if isempty(GVAR)
	 GVAR = [fieldnames(ptr) fieldnames(ptr)]' ;
      endif
      if PAR.udf
	 load(ufile = fullfile(ADIR, "udf", [GVAR{1,1} ".ot"])) ;
      else
	 Judf = false(columns(ptr.(GVAR{1,1}).x),1) ;
      endif
      p_ana.ptr = selfield(ptr, GVAR{1,:}) ;
      printf("plt: <-- %s\n", file = fullfile(CDIR, "ptrz.ob")) ;
      load(file) ; cur.ptr = selfield(ptr, GVAR{2,:}) ; ;
      printf("plt: <-- %s\n", file = fullfile(SDIR, "ptrz.ob")) ;
      load(file) ; fut.ptr = selfield(ptr, GVAR{2,:}) ;
      clear ptr ;
   endif

   for jv=1:length(VAR)

      v = VAR{jv} ; vn = VARNAME{jv} ; vv = strrep(v, ".", "_") ;
      if !any(J{jv}) || (exist("V", "var") && !strcmp(v, V)), continue ; endif
      vu = punit(v) ;
      JJ = find(strcmp(vars, v)')(J{jv}) ; udf = str2num(char(undefs(JJ)))' ;

      printf("plt: <-- %s\n", ifile = fullfile(ldata, [v ".csv"])) ;
      obs = read_stat(ifile) ;
      [JP JJ] = pick_meta(ids, vars, obs.ids, v) ;
      obs.ids = ids(JP)(JJ)' ;
      obs.vars = vars(JP)(JJ)' ;
      obs.names = names(JP)(JJ)' ;
      obs.lons = str2double(lons(JP)(JJ)') ;
      obs.lats = str2double(lats(JP)(JJ)') ;
      obs.alts = str2double(alts(JP)(JJ)') ;
      obs.udf = str2double(undefs(JP)(JJ)') ;
      obs.x = obs.x(:,JJ) ; obs.x = lim2nan(obs.x, obs.udf) ;
      TP.obs.(v).id = obs.id ; TP.obs.(v).x = obs.x ; TP.obs.(v).ll = [obs.lons ; obs.lats] ;
      ## Jv = cellfun(@(x) ~isempty(x), strfind(vars, v)) ;
      ## obs.names = tolower(names(Jv)) ; obs.vars = tolower(vars(Jv)) ; obs.cal = "gregorian" ;
      ## load(ofile = fullfile(ldata, "pdd.ob")) ;
      ## obs.lons = pdd.lons(J{jv}) ; obs.lats = pdd.lats(J{jv}) ;
      ## obs = pre_pare(pre_pare(obs), "inv") ;
 
      printf("plt: <-- %s\n", ifile = fullfile(ADIR, PAR.pdd, [vv ".csv"])) ;
      ana = read_stat(ifile) ;
      ana.x = lim2nan(ana.x, obs.udf) ;
      TP.ana.(v).id = ana.id ; TP.ana.(v).x = ana.x ; TP.ana.(v).ll = TP.obs.(v).ll ;

      printf("plt: <-- %s\n", ifile = fullfile(SDIR, PAR.pdd, [vv ".csv"])) ;
      sim1 = read_stat(ifile) ;
      sim1.x = lim2nan(sim1.x, obs.udf) ;
      TP.sim1.(v).id = sim1.id ; TP.sim1.(v).x = sim1.x ; TP.sim1.(v).ll = TP.obs.(v).ll ;

      printf("plt: <-- %s\n", ifile = fullfile(CDIR, PAR.pdd, [vv ".csv"])) ;
      sim2 = read_stat(ifile) ;
      sim2.x = lim2nan(sim2.x, obs.udf) ;
      TP.sim2.(v).id = sim2.id ; TP.sim2.(v).x = sim2.x ; TP.sim2.(v).ll = TP.obs.(v).ll ;

      if Lptr
	 [p_ana.Y wLptr] = ptr_loc(p_ana.ptr, vn, str2double(lons), str2double(lats)) ;
	 if wLptr
	    printf("plt: pairing %s with %s\n", vn, p_ana.Y.name) ;
	    printf("plt: using:\n") ;
	    printf(["\t" repmat("%d ", 1, length(p_ana.Y.K)) "\n"], p_ana.Y.K) ;
	    printf(["\t" repmat("%d ", 1, length(p_ana.Y.lon)) "\n"], p_ana.Y.lon) ;
	    printf(["\t" repmat("%d ", 1, length(p_ana.Y.lat)) "\n"], p_ana.Y.lat) ;
	    if (isT = strncmpi(vn, "temp", 4)) && mean(p_ana.Y.x) > 100, p_ana.Y.x = p_ana.Y.x - T0 ; endif
	    if isP = strncmpi(vn, "prec", 4)
	       if strcmp(PAR.ana, "ERAi")
		  p_ana.Y.x = p_ana.Y.x * 1000 ;
	       else
		  p_ana.Y.x = p_ana.Y.x * P0 ;
	       endif
	    endif
	    TP.p_ana.(v).id = p_ana.Y.id ; TP.p_ana.(v).x = p_ana.Y.x ; TP.p_ana.(v).ll = [p_ana.Y.lon' ; p_ana.Y.lat'] ;
	    cur.Y = ptr_loc(cur.ptr, vn, str2double(lons), str2double(lats)) ;
	    printf("plt: pairing %s with %s\n", vn, cur.Y.name) ;
	    if isT && mean(cur.Y.x) > 100, cur.Y.x = cur.Y.x - T0 ; endif
	    if isP, cur.Y.x = cur.Y.x * P0 ; endif
	    TP.cur.(v).id = cur.Y.id ; TP.cur.(v).x = cur.Y.x ; TP.cur.(v).ll = [cur.Y.lon' ; cur.Y.lat'] ;
	    fut.Y = ptr_loc(fut.ptr, vn, str2double(lons), str2double(lats)) ;
	    printf("plt: pairing %s with %s\n", vn, fut.Y.name) ;
	    if isT && mean(fut.Y.x) > 100, fut.Y.x = fut.Y.x - T0 ; endif
	    if isP, fut.Y.x = fut.Y.x * P0 ; endif
	    TP.fut.(v).id = fut.Y.id ; TP.fut.(v).x = fut.Y.x ; TP.fut.(v).ll = [fut.Y.lon' ; fut.Y.lat'] ;
	 endif
      endif
      ## continue ;

      for k=1:length(Sfun)
	 sfun = Sfun{k} ;
	 [obs.(sfun).id obs.(sfun).x] = annstat(obs.id, obs.x, sfun, annarg{:}) ;
	 [ana.(sfun).id ana.(sfun).x] = annstat(ana.id, ana.x, sfun, annarg{:}) ;
	 [sim1.(sfun).id sim1.(sfun).x] = annstat(sim1.id, sim1.x, sfun, annarg{:}) ;
	 [sim2.(sfun).id sim2.(sfun).x] = annstat(sim2.id, sim2.x, sfun, annarg{:}) ;
	 if wLptr
	    [p_ana.(sfun).id p_ana.(sfun).x] = annstat(p_ana.Y.id, p_ana.Y.x, sfun, annarg{:}) ;
	    [cur.(sfun).id cur.(sfun).x] = annstat(cur.Y.id, cur.Y.x, sfun, annarg{:}) ;
	    [fut.(sfun).id fut.(sfun).x] = annstat(fut.Y.id, fut.Y.x, sfun, annarg{:}) ;
	 endif
      endfor

      for j=1:columns(obs.x)
	 jj = JJ(j) ;
	 st = num2str(obs.ids{j}) ; 
	 st = strrep(st, " ", "_") ; wst = strrep(st, "_", "\\_") ;

	 ## annual timeseries
	 ##{
	 tsfile = fullfile(pdir, [v "." st ".ts.png"]) ;
	 tsfig = figure("visible", vis, "position", [0 1 0.2 0.2]) ;

	 for k=1:2
	    subplot(2, 1, k) ;
	    sfun = Sfun{k} ;
	    line(obs.(sfun).id(:,1), obs.(sfun).x(:,j), "linewidth", lw(1), "color", cmap(1,:)) ;
	    line(ana.(sfun).id(:,1), ana.(sfun).x(:,j), "linewidth", lw(2), "color", cmap(2,:)) ;
	    axis([1900 2100]) ;
	    line(sim2.(sfun).id(:,1), sim2.(sfun).x(:,j), "linewidth", lw(3), "color", cmap(3,:)) ;
	    line(sim1.(sfun).id(:,1), sim1.(sfun).x(:,j), "linewidth", lw(4), "color", cmap(4,:)) ;
	    lpdd = {"obs", PAR.ana, PAR.sim{2}, PAR.sim{1}} ;

	    if wLptr
	       p_ana.l = sprintf("%s (%.1f,%.1f)", PAR.ana, p_ana.Y.lon(jj), p_ana.Y.lat(jj)) ;
	       cur.l = sprintf("%s (%.1f,%.1f)", PAR.sim{2}, cur.Y.lon(jj), cur.Y.lat(jj)) ;
	       fut.l = sprintf("%s (%.1f,%.1f)", PAR.sim{1}, fut.Y.lon(jj), fut.Y.lat(jj)) ;
	       if str2num(strsplit(version, "."){2}) > 7
		  hl = legend(lpdd{:}, "location", "north") ;
		  legend("boxoff") ;
		  legend("right") ;
		  hptr = copyobj(gca, gcf) ; delete(get(hptr, "children")) ; delete(hl) ;
		  axes(hptr) ; set(hptr, "color", "none") ;
		  line(p_ana.(sfun).id(:,1), p_ana.(sfun).x(:,jj), "linewidth", lw(2)/2, "color", cmap(2,:), "linestyle", "--") ;
		  line(cur.(sfun).id(:,1), cur.(sfun).x(:,jj), "linewidth", lw(3)/2, "color", cmap(3,:), "linestyle", "--") ;
		  line(fut.(sfun).id(:,1), fut.(sfun).x(:,jj), "linewidth", lw(4)/2, "color", cmap(4,:), "linestyle", "--") ;
		  legend({p_ana.l, cur.l, fut.l}, "location", legpos) ;
		  legend("boxoff") ;
		  legend("right") ;
	       else
		  line(p_ana.(sfun).id(:,1), p_ana.(sfun).x(:,jj), "linewidth", lw(2)/2, "color", cmap(2,:), "linestyle", "--") ;
		  line(cur.(sfun).id(:,1), cur.(sfun).x(:,jj), "linewidth", lw(3)/2, "color", cmap(3,:), "linestyle", "--") ;
		  line(fut.(sfun).id(:,1), fut.(sfun).x(:,jj), "linewidth", lw(4)/2, "color", cmap(4,:), "linestyle", "--") ;
		  [hleg, hleg_obj, hplot, labels] = legend({lpdd{:}, p_ana.l, cur.l, fut.l}, "location", legpos) ;
		  legend("boxoff") ;
		  legend("right") ;
	       endif
	    else
	       legend(lpdd, "location", legpos) ;
	       legend("boxoff") ;
	       legend("right") ;
	    endif
	    xlabel("year") ; ylabel(["annual " sfun " " v vu]) ;
	    title([upper(PAR.pdd(1)) PAR.pdd(2:end) ", " v ", " wst]) ;
	 endfor

	 printf("plt: --> %s\n", tsfile) ;
	 print(tsfile, "-r600") ;
	 ##pause(10) ;
	 close(tsfig) ;
	 ##}

	 ##{
	 ## boxplots
	 bpfile = fullfile(pdir, [v "." st ".bp.png"]) ;
	 bpfig = figure("visible", vis, "position", [0 1 0.2 0.2]) ;
	 box.xm = {} ;
	 for m=1:12
	    obs.I = obs.id(:,2) == m ;
	    ana.I = ana.id(:,2) == m ;
	    sim1.I = sim1.id(:,2) == m ;
	    box.xm{1}{m} = obs.x(obs.I,j) ;
	    box.xm{2}{m} = ana.x(ana.I,j) ;
	    sim2.I = sim2.id(:,2) == m ;
	    box.xm{3}{m} = sim2.x(sim2.I,j) ;
	    box.xm{4}{m} = sim1.x(sim1.I,j) ;
	 endfor
	 hg = aboxplot(box.xm, "labels", {mth{:}}, "colormap", cmap) ;
	 legend(hg, {"obs", PAR.ana, PAR.sim{2}, PAR.sim{1}}, "location", legpos) ;
	 legend("boxoff") ;
	 legend("right") ;
	 xlabel("month") ; ylabel(v) ;
	 title([upper(PAR.pdd(1)) PAR.pdd(2:end) ", " v ", " wst]) ;
	 printf("plt: --> %s\n", bpfile) ;
	 print(bpfile, "-r600", "-portrait") ;
	 close(bpfig) ;
	 ##}

	 ##{
	 dfile = fullfile(pdir, [v "." st ".d.png"]) ;
	 #### workaround for fltk
	 ##gt = graphics_toolkit ; graphics_toolkit fltk ;
	 ##dfig = figure("visible", "on") ;
	 ####pause(10) ;
	 ##close(dfig) ;
	 ## end workaround for fltk
	 dfig = figure("visible", vis, "position", [0 1 0.2 0.2]) ;
	 wlw = lw ; if BATCH, wlw = lw / 5 ; endif
	 plt_dly(dfig, obs, ana, j, v, vn, per, wlw, cmap) ;
	 printf("plt: --> %s\n", dfile) ;
	 print(dfile, "-r600", "-portrait") ;
	 ##pause(10) ;
	 close(dfig) ;
	 ##graphics_toolkit(gt) ;
	 ##}

      endfor

   endfor

   if !isfield(TP.obs, "T") || !isfield(TP.obs, "P"), continue ; endif

   ## T-P kernel regression
   nb = 10 ;
   [obs.T obs.Pf obs.Pi] = TP_kreg(TP.obs, nb) ;
   [ana.T ana.Pf ana.Pi] = TP_kreg(TP.ana, nb) ;
   [sim2.T sim2.Pf sim2.Pi] = TP_kreg(TP.sim2, nb) ;
   [sim1.T sim1.Pf sim1.Pi] = TP_kreg(TP.sim1, nb) ;
   if wLptr
      [p_ana.T p_ana.Pf p_ana.Pi] = TP_kreg(TP.p_ana, nb) ;
      [cur.T cur.Pf cur.Pi] = TP_kreg(TP.cur, nb) ;
      [fut.T fut.Pf fut.Pi] = TP_kreg(TP.fut, nb) ;
   endif

   for j=1:columns(obs.Pf)

      jj = JJ(j) ;
      st = num2str(ids(jj){:}) ; 
      st = strrep(st, " ", "_") ; wst = strrep(st, "_", "\\_") ;

      TPfile = fullfile(pdir, [v "." st ".TP.png"]) ;
      TPfig = figure("visible", vis) ;

      subaxis(2, 1, 1) ;
      line(obs.T{j}, obs.Pf{j}, "linewidth", lw(1), "color", cmap(1,:)) ;
      line(ana.T{j}, ana.Pf{j}, "linewidth", lw(2), "color", cmap(2,:)) ;
      line(sim2.T{j}, sim2.Pf{j}, "linewidth", lw(3), "color", cmap(3,:)) ;
      line(sim1.T{j}, sim1.Pf{j}, "linewidth", lw(4), "color", cmap(4,:)) ;
      title([upper(PAR.pdd(1)) PAR.pdd(2:end) ", " v ", " wst]) ;
      if wLptr
	 line(p_ana.T{j}, p_ana.Pf{j}, "linewidth", lw(2)/2, "color", cmap(2,:), "linestyle", "--") ;
	 line(cur.T{j}, cur.Pf{j}, "linewidth", lw(2)/2, "color", cmap(3,:), "linestyle", "--") ;
	 line(fut.T{j}, fut.Pf{j}, "linewidth", lw(2)/2, "color", cmap(4,:), "linestyle", "--") ;
	 [hleg, hleg_obj, hplot, labels] = legend({lpdd{:}, p_ana.l, cur.l, fut.l}, "location", legpos) ;
	 legend("boxoff") ;
	 legend("right") ;
      else
	 legend(lpdd, "location", legpos) ;
	 legend("boxoff") ;
	 legend("right") ;
      endif
      ylabel(["Pf"]) ;

      subaxis(2, 1, 2) ;
      line(obs.T{j}, obs.Pi{j}, "linewidth", lw(1), "color", cmap(1,:)) ;
      line(ana.T{j}, ana.Pi{j}, "linewidth", lw(2), "color", cmap(2,:)) ;
      line(sim2.T{j}, sim2.Pi{j}, "linewidth", lw(3), "color", cmap(3,:)) ;
      line(sim1.T{j}, sim1.Pi{j}, "linewidth", lw(4), "color", cmap(4,:)) ;
      if wLptr
	 line(p_ana.T{j}, p_ana.Pi{j}, "linewidth", lw(2)/2, "color", cmap(2,:), "linestyle", "--") ;
	 line(cur.T{j}, cur.Pi{j}, "linewidth", lw(2)/2, "color", cmap(3,:), "linestyle", "--") ;
	 line(fut.T{j}, fut.Pi{j}, "linewidth", lw(2)/2, "color", cmap(4,:), "linestyle", "--") ;
	 [hleg, hleg_obj, hplot, labels] = legend({lpdd{:}, p_ana.l, cur.l, fut.l}, "location", legpos) ;
	 legend("boxoff") ;
	 legend("right") ;
      else
	 legend(lpdd, "location", legpos) ;
	 legend("boxoff") ;
	 legend("right") ;
      endif
      xlabel("T [C]") ;
      ylabel(["Pi [mm/d]"]) ;
      ## orient landscape
      ## set(gca, "xgrid", "on") ;
      set(gcf, "paperposition", [0.05 0.05 0.95 0.95])
      printf("plt: --> %s\n", TPfile) ;
      print(TPfile, "-r600") ;
      close(TPfig) ;

   endfor

endfunction


function [y rv] = ptr_loc (ptr, vn, lon, lat)
 
   ## usage: [y rv] = ptr_loc (ptr, vn, lon, lat)
   ##
   ## nearest gridpoint for vn
   ## UGLY

   global ADIR GVAR Judf

   y = [] ;
   vn = tolower(vn) ;
   if !isempty(strfind(vn, "temp")), vn = "temp" ; endif
   if !isempty(strfind(vn, "precip")), vn = "precip" ; endif

   for [v k] = ptr
      name = tolower(v.name) ;
      if !isempty(strfind(name, "temp")), name = "temp" ; endif
      if !isempty(strfind(name, "wind")), name = "wind" ; endif
      if !isempty(strfind(name, "sun")), name = "sun" ; endif
      if !isempty(strfind(name, "precip")), name = "precip" ; endif
      if strcmp(name, vn)
	 clear y ;
	 y.name = name ;
	 break ;
      endif
   endfor

   if !(rv = isfield(y, "name"))
      return ;
   endif

   if isfield(v, "lon")
      xlon = v.lon ;
   else
      xlon = lon ;
   endif
   if isfield(v, "lat")
      xlat = v.lat ;
   else
      xlat = lat ;
   endif

   if isvector(xlon) && isvector(xlat)
      [xlon xlat] = meshgrid(xlon, xlat) ;
   endif

   v.x = v.x(:,!Judf) ; xlon = xlon(!Judf) ; xlat = xlat(!Judf) ;

   y.id = v.id ; y.name = k ;
   y.K = findgp([lon(:)' ; lat(:)'], [xlon(:)' ; xlat(:)']) ;
   y.lon = xlon(:)(y.K) ; y.lat = xlat(:)(y.K) ; y.x = v.x(:,y.K) ;

endfunction


function v = selfield (u, varargin)
	 
   ## usage: v = selfield (u, varargin)
   ## 
   ## select fields from structure

   for [f k] = u
      for i = 2:nargin
	 if !strcmp(varargin{i-1}, k), continue ; endif
	 v.(k) = f ;
      endfor
   endfor

endfunction


function plt_dly (h, obs, ana, j, v, vn, per, lw, cmap)

   ## usage: plt_dly (h, obs, ana, j, v, vn, per, lw, cmap)
   ## 
   ## plot daily values

   global BATCH PAR

   vis = get(h, "Visible") ;
   figure(h) ;

   st = num2str(obs.ids{j}) ; 
   st = strrep(st, " ", "_") ; wst = strrep(st, "_", "\\_") ;
   vu = punit(v) ;

   obs.I = find(sdate(obs.id, per)) ;
   ana.I = find(sdate(ana.id, per)) ;

   if L = strcmp(graphics_toolkit, "fltk")
      dfun = @dateax ;
   else
      dfun = @pltid ;
   endif

   oI = obs.I(1:floor(length(obs.I)/2)) ; aI = ana.I(1:floor(length(ana.I)/2)) ;
   line(dfun(obs.id(oI,:)), obs.x(oI,j), "color", cmap(1,:), "linewidth", lw(1)) ;
   ## line(dfun(ana.id(aI,:)), ana.x(aI,j), "color", cmap(2,:), "linewidth", lw(2)) ;
   oI = obs.I(ceil(length(obs.I)/2):end) ; aI = ana.I(ceil(length(ana.I)/2):end) ;
   line(dfun(obs.id(oI,:)), obs.x(oI,j), "color", cmap(1,:), "linewidth", lw(1)) ;
   ## line(dfun(ana.id(aI,:)), ana.x(aI,j), "color", cmap(2,:), "linewidth", lw(2)) ;
   xlim = get(gca, "xlim") ;
   ax = axis ; clf ; set(gcf, "Visible", vis) ;

   ##subaxis(2, 1, 1) ;
   subplot(2, 1, 1) ;
   oI = obs.I(1:floor(length(obs.I)/2)) ; aI = ana.I(1:floor(length(ana.I)/2)) ;
   if any(~isnan(obs.x(oI,j))), axis(ax) ; endif
   line(dfun(obs.id(oI,:)), obs.x(oI,j), "color", cmap(1,:), "linewidth", lw(1)) ;
   line(dfun(ana.id(aI,:)), ana.x(aI,j), "color", cmap(2,:), "linewidth", lw(2)) ;
   if L
      datetick("dd mmm yy") ;
   else
      set(gca, "xlim", [xlim(1) mean(xlim)]) ;
   endif
   legend({"obs" "ERAi"}, "location", "northwest") ;
   legend("boxoff") ; legend("right") ;
   title([upper(PAR.pdd(1)) PAR.pdd(2:end) ", " vn ", " wst]) ;
   ylabel([v vu]) ;

   ##subaxis(2, 1, 2) ;
   subplot(2, 1, 2) ;
   oI = obs.I(ceil(length(obs.I)/2):end) ; aI = ana.I(ceil(length(ana.I)/2):end) ;
   if any(~isnan(obs.x(oI,j))), axis(ax) ; endif
   line(dfun(obs.id(oI,:)), obs.x(oI,j), "color", cmap(1,:), "linewidth", lw(1)) ;
   line(dfun(ana.id(aI,:)), ana.x(aI,j), "color", cmap(2,:), "linewidth", lw(2)) ;
   if L
      datetick("dd mmm yy") ;
   else
      set(gca, "xlim", [mean(xlim) xlim(2)]) ;
   endif
   legend({"obs" "ERAi"}, "location", "northwest") ;
   legend("boxoff") ; legend("right") ;
   ylabel([v vu]) ;

endfunction


function vu = punit (v)

   ## usage: vu = punit (v)
   ## 
   ## physical units

   switch v
      case "T"
	 vu = " [C]" ;
      case "P"
	 vu = " [mm/d]" ;
      otherwise
	 vu = "" ;
   endswitch

endfunction


function [Tq Pf Pi] = TP_kreg (s, nb)

   ## usage: [Tq Pf Pi] = TP_kreg (s, nb)
   ## 
   ## T-P kernel regression

   binn = "T" ;

   [IT IP] = common(s.T.id, s.P.id) ;
   s.T.id = s.T.id(IT,:) ; s.T.x = s.T.x(IT,:) ;
   s.P.id = s.P.id(IP,:) ; s.P.x = s.P.x(IP,:) ;

   ## I = ssn(s.T.id, "cold") ;
   ## I = ssn(s.T.id, "warm") ;
   I = ssn(s.T.id, "all") ;

   for j = 1:columns(s.P.ll)

      s.P.T(j) = findll(s.P.ll(:,j), s.T.ll) ;
      X.id = s.P.id(I,:) ;
      X.x = [s.T.x(I,s.P.T(j)) s.P.x(I,j)] ;

      if 0
	 scatter(X.x(:,1), X.x(:,2), 5, "black", "filled") ;
	 xlabel("T [C]") ; ylabel(["P [mm]"])
	 orient landscape
	 set(gca, "xgrid", "on") ;
	 set(gcf, "paperposition", [0.05 0.05 0.95 0.95])
	 print(fullfile(cloud, ["T-P_" sfx "_scatter.png"])) ;
      endif

      ## pre-select data
      II = all(!isnan(X.x), 2) ;
      x = X.x(II,:) ;
      Iw = x(:,2) > 0 ;	# wet days

      ## quantiles
      w = quantile(x(:,1), linspace(0, 1, nb+1))' ;
      Tq{j} = mean([w(1:end-1)', w(2:end)'], 2) ;
      ## Tq{j} = linspace(min(x(:,1)), max(x(:,1)), nb+1)' ;

      ## frequencies
      IT = arrayfun(@(u,v) u < x(:,1) & x(:,1) < v, w(1:end-1), w(2:end), "UniformOutput", false) ;
      Pf{j} = cell2mat(cellfun(@(i) sum(Iw&i)/sum(i), IT, "UniformOutput", false))' ;

      ## intensity (kernel regression)
      n = rows(x) ;
      Pi{j} = kernel_regression(Tq{j}, x(Iw,2), x(Iw,1), n^(-1/5), "kernel_normal", true, false) ;
      ## Pi{j} = qntl(x(Iw,:), binn, 0.5, nb+1) ;

   endfor

endfunction
